Double Negative Control Inference in Test-Negative Design Studies of Vaccine Effectiveness

The test-negative design (TND) has become a standard approach to evaluate vaccine effectiveness against the risk of acquiring infectious diseases in real-world settings, such as Influenza, Rotavirus, Dengue fever, and more recently COVID-19. In a TND study, individuals who experience symptoms and seek care are recruited and tested for the infectious disease which defines cases and controls. Despite TND's potential to reduce unobserved differences in healthcare seeking behavior (HSB) between vaccinated and unvaccinated subjects, it remains subject to various potential biases. First, residual confounding may remain due to unobserved HSB, occupation as healthcare worker, or previous infection history. Second, because selection into the TND sample is a common consequence of infection and HSB, collider stratification bias may exist when conditioning the analysis on tested samples, which further induces confounding by latent HSB. In this paper, we present a novel approach to identify and estimate vaccine effectiveness in the target population by carefully leveraging a pair of negative control exposure and outcome variables to account for potential hidden bias in TND studies. We illustrate our proposed method with extensive simulations and an application to study COVID-19 vaccine effectiveness using data from the University of Michigan Health System.


Contents
1 Introduction

Text-negative design studies of vaccine effectiveness
The test-negative design (TND) has become a standard approach to evaluate real-world vaccine effectiveness (VE) against the risk of acquiring infections diseases (Chung et al., 2020;Flannery et al., 2019;Jackson et al., 2017;Rolfes et al., 2019;Tenforde et al., 2021). In an outpatient Influenza VE test-negative design, for example, symptomatic individuals seeking care and meeting eligibility criteria are enrolled and their Influenza virus infection status is subsequently confirmed via a laboratory test. VE against flu infection is then measured by comparing the prevalence of vaccination between the test-positive "cases" and test-negative "controls" (Jackson and Nelson, 2013;Jackson et al., 2017). Besides Influenza, the TND and its variants have also been applied to study VE against pneumococcal disease (Broome, Facklam, and Fraser, 1980), dengue (Anders et al., 2018;Utarini et al., 2021), rotavirus (Boom et al., 2010;Schwartz et al., 2017), and other infectious diseases. Recently, the TND has increasingly been used in post-licensure evaluation of COVID-19 VE (Dagan et al., 2021;Dean, Hogan, and Schnitzer, 2021;Hitchings et al., 2021;Olson et al., 2022;Patel, Jackson, and Ferdinands, 2020;Thompson et al., 2021).
Test-negative designs are believed to reduce unmeasured confounding bias due to healthcare seeking behavior (HSB), whereby care seekers are more likely to be vaccinated, have healthier behaviors that reduces the risk of infection, and get tested when ill (Jackson et al., 2006;Shrank, Patrick, and Brookhart, 2011). By restricting analysis to care seekers who are tested for the infection in view (e.g. Influenza or , the vaccinated and unvaccinated are more likely to share similar HSB and underlying health characteristics. Misclassification of infection status is also reduced because the analysis is restricted to tested individuals (Jackson and Nelson, 2013). Sullivan et al. (2016) used directed acyclic graphs (DAG) to illustrate the rationale behind TND in the context of evaluating VE against Influenza infection, as shown in Figures 1(a) and (b). We denote Influenza vaccination status by A and Influenza infection by Y , so that the arrow A → Y represents VE against flu infection. Selection into the TND study sample, denoted by S, is triggered by a subject experiencing flu-like symptoms or acute respiratory illness, seeking care at clinics or hospitals, and getting tested for Influenza infection, hence the Y → S edge.
Healthcare seeking behavior, denoted by HSB, may affect S, A, and Y because subjects with certain health seeking proclivities may be more likely to seek care, take annual flu shots, and participate in healthy and preventative behaviors. The above variables are subject to effects of other clinical or demographic factors, such as age, season and high-risk conditions, included in Figure 1 as confounders X. The TND assumes that by restricting recruitment to care seekers, the study subjects have identical healthcare seeking behavior; in other words, conditioning the analysis on S = 1 necessarily leads to HSB= 1, which blocks the effects of HSB (Figure 1(b)).
The effects of X are further adjusted for by including these factors in a logistic regression model or by inverse probability weighting (Bond, Sullivan, and Cowling, 2016;Thompson et al., 2021).
However, the TND remains subject to potential hidden bias. First, the assumption that all study subjects seeking care are lumped into a single category HSB= 1 may be unrealistic. It may be more realistic that HSB is not a deterministic function of S and remains a source of confounding bias even after conditioning on S. Furthermore, there might be other mismeasured or unmeasured confounders, denoted as U . For example, healthcare workers are at increased risk of flu infection due to higher exposure to flu patients and are more likely to seek care and receive vaccination due to health agency guidelines (Black et al., 2018). Previous flu infection history may also be a source of confounding if it alters the likelihood of vaccination and care seeking, while also providing immunity against circulating strains (Krammer, 2019;Sullivan, Tchetgen Tchetgen, and Cowling, 2016). These unmeasured or mismeasured potential sources of confounding, if not properly accounted for, can result in additional confounding bias, as illustrated in Figure 1(c). Finally, collider stratification bias is likely present due to conditioning on S, which is a common consequence of HSB, other risk factors (X, U ), and Influenza infection Y (Lipsitch, Jha, and Simonsen, 2016). That is, conditioning on S unblocks the backdoor path A ← (X, U, HSB) → S ← Y , which would in principle be blocked if study subjects had identical levels of HSB and other risk factors (Sullivan, Tchetgen Tchetgen, and Cowling, 2016).
Accounting for these potential sources of bias is well known to be challenging, and potentially infeasible without additional assumptions or data. This can be seen in Figure Figure 1: Causal relationships of variables in a test-negative design. Sullivan, Tchetgen Tchetgen, and Cowling, 2016 used (a) to illustrate the causal relationship between variables in a test-negative design in the general population, and used (b) to illustrate the assumption implicit in the common approach to estimate VE from the study data that study subjects have identical healthcare seeking behavior (HSB) (Sullivan, Tchetgen Tchetgen, and Cowling, 2016
The framework requires that at least one of two types of negative control variables are available which are a priori known to satisfy certain conditions: a negative control exposure (NCE) known to have no direct effect on the primary outcome; or a negative control outcome (NCO), known not to be an effect of the primary exposure. Such negative control variables are only valid and therefore useful to address unmeasured confounding in a given setting to the extent that they are subject to the same source of confounding as the exposure-outcome relationship of primary interest. Thus, the observed association between a valid NCE and the primary outcome (conditional on the primary treatment and observed covariates) or that between a valid NCO and the primary exposure can indicate the presence of residual confounding bias. For example, in a cohort study to investigate flu VE against hospitalization and death among seniors, to detect the presence of confounding bias due to underlying health characteristics, Jackson et al. (2006) used hospitalization/death before and after the flu season as NCOs and found that the association between flu vaccination and hospitalization was virtually the same before and during the flu season, suggesting that the lower hospitalization rate observed among vaccinated seniors versus unvaccinated seniors was partially due to healthy user bias.
Recently, new causal methods have been developed to not only detect residual confounding when present, but also to potentially de-bias an observational estimate of a treatment causal effect in the presence of unmeasured confounders when both an NCE and an NCO are available, referred to as the double negative control (Miao, Geng, and Tchetgen Tchetgen, 2018;Shi et al., 2020;Tchetgen Tchetgen et al., 2020). In this recent body of work, the double negative control design was extended in several important directions including settings in which proxies of treat-ment and outcome confounding routinely measured in well designed observational studies may be used as negative control variables, a framework termed proximal causal inference; longitudinal settings where one is interested in the joint effects of time-varying exposures (Ying et al., 2021), potentially subject to both measured and unmeasured confounding by time-varying factors; and in settings where one aims to estimate direct and indirect effects in mediation analysis subject to unmeasured confounding or unmeasured mediators (Dukes, Shpitser, and Tchetgen Tchetgen, 2021;Ghassami, Shpitser, and Tchetgen Tchetgen, 2021). Additional recent papers in this fast-growing literature include Qi, Miao, and Zhang (2021), Liu and Tchetgen Tchetgen (2021), Egami and Tchetgen Tchetgen (2021), Kallus, Mao, and Uehara (2021), Imbens, Kallus, and Mao (2021), Deaner (2018), Deaner (2021), Ghassami et al. (2021), Mastouri et al. (2021) and Ghassami, Shpitser, and Tchetgen Tchetgen (2022). Importantly, existing identification results in negative control and proximal causal inference literature has been restricted to i.i.d settings (Miao, Shi, and Tchetgen Tchetgen, 2018) and time series settings (Shi et al., 2021), and to date, to the best of our knowledge, outcome-dependent sampling settings such as TND, have not been considered, particularly one where confounding and selection bias might co-exist.

Outline
The rest of the paper is organized as followed: we introduce notation and the identification challenge in view in Sections 2.1. Next we develop the identification strategy and describe a new debiased estimator under a double negative control TND study in Section 2.2-2.4, assuming (1) homogeneous VE across strata defined by all measured and unmeasured confounders and (2) no direct effect of vaccination on selection into the TND sample. In Section 2.5, we relax the homogeneous VE assumption and describe identification and estimation allowing for VE to depend on observed covariates. In Section 2.6, we relax the assumption of no direct effect of vaccination on selection and introduce the assumptions under which our VE estimator is unbiased on the odds ratio scale. In Section 3, we demonstrate the performance of our method with simulation. In Section 4, the approach is further illustrated in an application to estimate COVID-19 VE against infection in a TND study nested within electronic health records from University of Michigan Health System. We then conclude with a discussion in Section 5. To fix ideas, we first review estimation assuming all confounders (U, X) are fully observed and the study sample is randomly drawn (rather than selected by testing) from source population, referred to as the "target population". That is, we observe data on (A, Y, U, X) which are independent and identically distributed in the target population. For each individual, we write Y (a) as the binary potential infection outcome had, possibly contrary to fact, the person's vaccination status been A = a, a = 0, 1. Our goal is to provide identification and estimation strategies for the causal risk ratio (RR) defined as Let β 0 denote the log causal RR, i.e., RR = exp(β 0 ). Following Hudgens and Halloran (2006) and Struchiner and Halloran (2007), we define VE as one minus the causal RR: V E = 1 − RR = 1 − exp(β 0 ). The potential outcomes and the observed data are related through the following assumptions: Assumption 1. (Identification conditions of mean potential outcomes).
Assumption 1(a) states that the infection status of a subject with vaccination status A = a is equal to the corresponding potential outcome Y (a). This further requires that the treatment is sufficiently well-defined and a subject's potential outcome is not affected by the treatment of other subjects (Cole and Frangakis, 2009). Assumption 1(b) states that treatment is exchangeable within strata of (U, X), i.e. there is no unmeasured confounding given (U, X). We develop methods that allow U to be unmeasured in Section 2.3. Assumption 1(c) states that for all realized values of (U, X) there is at least one individual with an opportunity to get treatment a = 0, 1.
Let Q(A = a, U, X) = 1/P (A = a|U, X) denote the inverse of the probability of vaccination status A = a given confounders (Rosenbaum, 1987). Under Assumption 1, it is well known that, if U were observed, the mean potential outcome for treatment a in the general population can be identified by inverse probability of treatment weighting (IPTW): for a = 0, 1. Therefore, the log causal RR β 0 satisfies the following equation Equivalently, we have for the unbiased estimating function, where V 0 (A, Y, U, X; β) = (−1) 1−A Q(A, U, X)Y exp(−βA).

Tackling selection bias under a semiparametric risk model
Next, consider a TND study for which data (A, Y, X, U ) is observed only for the tested individuals with S = 1. Because S is impacted by other factors such as infection, the estimating function V 0 (A, Y, U, X; β 0 ) may not be unbiased with respect to the study sample; i.e.
without another assumption about the selection process into the TND sample.
For a study sample of size n from a TND, we denote the i-th study subject's variables . . , n. For generalizability, we first make the key assumption that vaccination A is unrelated to selection S other than through a subject's infection status Y and confounders (U, X).
In the test-negative design, this assumption requires that an individual's decision to seek care and get tested only depends on the presence of symptoms and his/her underlying behavioral or socioeconomic characteristics, including HSB (contained in (U, X)); a person's vaccination status does not directly Influence their selection process. The DAGs in Figure 1(a)-(e) in fact encode this conditional independence condition. We will relax this assumption in Section 2.6.
Assumption 3 defines a semiparametric multiplicative risk model which posits that vaccine effectiveness, measured on the RR scale, is constant across (U, X) strata in the target population.
In other words, the effect of vaccination A on the risk of infection Y is not modified by confounders U, X. In Section 2.5, we will relax the assumption to allow for effect modification by measured The proof of Proposition 1 is in Appendix A. From Proposition 1, the IPTW estimating function V 0 derived from the target population is also unbiased with respect to the study sample.
Under Assumptions 1-3, one can estimate β 0 withβ as the solution to where n is the size of the selected sample, is the estimated probability of have vaccination status A = A i given confounders (U i , X i ). The resulting estimator is essentially the IPTW estimator of marginal RR in Schnitzer (2022) assuming (U i , X i )'s are all observed.
However, Q(A, U, X) cannot be estimated because U is unobserved. Furthermore, even if U were observed, Q(A i , U i , X i ) may not be identified from the TND sample due to selection bias.
In the next section, we describe a new framework to account for unmeasured confounding in a TND setting, leveraging negative control exposure and outcome variables.
2.3 Tackling unmeasured confounding bias leveraging negative controls 2.3.1 Negative control exposure (NCE) and treatment confounding bridge function As shown in Figure 1(e), suppose that one has observed a valid possibly vector-valued NCE, denoted as Z, which is a priori known to satisfy the following key independence conditions: Assumption 4 essentially states that any existing Z − Y association conditional on (X, A) in the target population must be a consequence of their respective association with U , therefore indicating the presence of confounding bias. Importantly, the NCE must a priori be known to have no causal effect on infection status (Miao, Shi, and Tchetgen Tchetgen, 2018). Likewise, the association between Z and S conditional on (X, A) is completely due to their respective association with U . Figure 1(e) presents a graphical illustration of an NCE that satisfies Assumption 4.
In the Influenza VE setting, a candidate NCE can be vaccination status for the preceding year, or other vaccination status such as Tdap (Tetanus, Diphtheria, Pertussis) vaccine, as both are known to effectively provide no protection against the circulating flu strain in a given year.
We now provide an intuitive description of our approach to leverage Z as an imperfect proxy of U for identification despite being unable to directly observe U .
To illustrate the rationale behind identification, ignore selection bias for now and suppose that Q(A, U ) = α 0 + α 1 A + α 2 U , suppressing measured confounders X. Although U is unobserved, suppose further that Z satisfies E[Z|A, U ] = γ 0 + γ 1 A + γ 2 U . Then we have that would naturally follow that the IPTW method in (1) can be recovered by Therefore, β 0 can be identified if the distribution of (A, Y, Z) in the target population is available provided that parameters indexing q can be identified.
The above insight motivates the following assumption: Assumption 5 (treatment confounding bridge function). There exists a function q(A, Z, X) that satisfies, for every a, u and x, The function q that satisfies (7) is called a treatment confounding bridge function, as it bridges the observed NCE with the unobserved propensity score (Cui et al., 2020). Below we give two examples where the integral equation (7) can easily be solved and the treatment confounding bridge function q admits a closed form solution.
Example 1. (Binary U and Z) Suppose that U is binary, and so is the NCE Z. For simplicity we suppress X. The integral equation (7) can then be written as 1 z=0 q(a, z)P (Z = z|U = u, A = a) = P (A = a|U = u) −1 , or equivalently, 1 z=0 p za.u q(a, z) = 1 for each a, u ∈ {0, 1}, where p za.u = P (Z = z, A = a|U = u). Therefore, the treatment confounding bridge function q(a, z) solves the linear equation If the matrix P Z,A|U is invertible, then q(a, z) has a closed form solution given by The result can be extended to the cases where Z is polytomous as detailed in Appendix D.
Example 2. (Continuous U and Z) Suppose the unmeasured confounder U and the NCE Z are continuous. Further assume that By the derivation in Appendix E, the treatment confounding bridge function q(A, Z, X) is Formally, Equation (7) defines a Fredholm integral equation of the first kind, with treatment confounding bridge function q(A, Z, X) as its solution (Cui et al., 2020). Heuristically, the existence of a treatment confounding bridge function requires that variation in Z induced by U is sufficiently correlated with variation in A induced by U . For instance, in Example 1, existence of a treatment confounding bridge function requires that the matrix P Z,A|U is nonsingular. In Example 2, the existence of a treatment confounding bridge function amounts to the condition µ U Z = 0, which again requires Z ⊥ ⊥ U |A, X. Cui et al. (2020) provided formal conditions sufficient for the existence of the treatment confounding bridge function satisfying Equation (7). These conditions are reproduced for completeness in Appendix B.
Thus, under Assumption 5, we propose to construct a new unbiased estimating function for Theorem 1. (Moment restriction of β 0 ) Under Assumptions 1-5, we have that The proof of Theorem 1 is in Appendix C. In practice, if one can consistently estimate the treatment confounding bridge function q(A, Z, X) with q(A, Z, X), Theorem 1 suggests estimating β 0 by solving the estimating equation which results in a closed form estimator Importantly, although (7) may not have a unique solution, any solution uniquely identifies the causal log RR β 0 . The result in Theorem 1 cannot directly be applied in practice because the treatment confounding bridge function is not identifiable even if random samples from the target population were available -solving (7) requires additional information about U which is unobserved. For instance, in Example 1 one is unable to directly estimate q(a, z) because p za.u in (8) cannot directly be estimated from the observed data.

Negative control outcome (NCO) for identification of treatment confounding bridge function
For identification and estimation of q, we leverage negative control outcomes (NCO) to construct feasible estimating equations for the treatment confounding bridge function as in Cui et al. (2020). Similar to NCEs, NCOs can be viewed as imperfect proxies of U . However, unlike NCEs, a valid NCO, denoted by W , is a measured covariate which is (i) known a priori not to be a causal effect of either the primary exposure A nor the NCE Z; and (ii) is associated with (A, Z) conditional on X only to the extent that it is associated with U .
Formally, we make the following assumption.
Similar to Cui et al. (2020), we leverage the availability of an NCO as an additional proxy to identify the treatment confounding bridge function. However, a complication arises due to lack of a random sample from the target population, a key requirement in the approach outlined in Cui et al. (2020). In general, it is not possible to obtain sufficient information about neither the distribution of W nor that of U in the target population from the TND data without an additional structural assumption (Bareinboim and Pearl, 2012). In the following, we avoid imposing such an additional structural assumption by leveraging an important feature of infectious diseases such as Influenza and COVID-19; mainly that contracting such an infection is a rare event in most target populations of interest, and therefore information from the target population that is relevant for estimating the treatment confounding bridge function can be recovered from the test-negative control group. Formally, we make the following rare disease assumption.

Assumption 7 (Rare infection).
There exist a small positive number δ > 0 such that Assumption 7 states that infected subjects, whether vaccinated or not and regardless of their negative control outcomes, only constitute a small proportion of each (U, X) stratum in the general population; specifically, the assumption implies that 1 1−δ ≤ P (A,Z|U =u,X=x,Y =0) Thus, under Assumptions 2, 4 and 7, P (A = a, Z = z|U = u, X = x) ≈ P (A = a, Z = z|U = u, X = x, Y = 0, S = 1) for all a, z, x, u. We now introduce a key property of the treatment confounding bridge function in Theorem 2, which is proved in Appendix F.
Theorem 2 (Identification of the treatment confounding bridge function). Under Assumptions 2, 4, 5, 6, and 7, for a = 0, 1 we have that Thus, provided δ ≈ 0, Theorem 2 suggests that an approximation to the treatment confounding bridge function can be obtained by solving the following integral equation involving only observed data as long as a solution exists. Accordingly, hereafter suppose that the following assumption holds.
Assumption 8 (Existence of a unique solution to (12)). There exists a unique square-integrable function q * (A, Z, X) that satisfies (12).
Heuristically, uniqueness of a solution to (12) requires that variation in W is sufficiently informative about variation in Z, in the sense that there is no source of variation in W that is not associated with a corresponding source of variation in Z. See Appendix G for further elaboration of completeness conditions and D'Haultfoeuille (2011) and Newey and Powell (2003) for related use of the assumption in the literature. Below we briefly illustrate Assumption 8 in the examples of Section 2.3.1.
Example 1 . Suppose U and Z are both binary, and a binary NCO W is also observed. Let (12) is equivalent to solving the system of linear equations p 0a.0 q * (a, 0) + p 1a.0 q * (a, 1) = 1; p 0a.1 q * (a, 0) + p 1a.1 q * (a, 1) = 1, that the probabilities p za.w can all be estimated from the study sample.
We emphasize that the solution to Equation (12) is ultimately an approximation to the (non-identifiable) treatment confounding bridge function in the target population. The accuracy of this approximation relies on the extent to which the rare disease assumption holds in the target population of interest. We study the potential bias resulting from a departure of this key assumption in the Appendix I. We further observe that, under the null hypothesis of no vaccine effectiveness, or if W has no direct effects on Y or S, then the function q * (A, Z, X) equals the treatment confounding bridge exactly, even for a non-rare disease outcome, as stated in the corollary below. We prove Corollary 1 in Appendix F. From Theorem 2, we immediately have the following corollary which provides a basis for estimation of q * (A, Z, X) from the observed TND sample.
Corollary 2. Under the conditions listed in Theorem 2, for any function m(W, A, X), the solution q * (A, Z, X) to Equation (12) also solves the population moment equation We prove Corollary 2 in Appendix H. In practical situations where a parametric model q * (A, Z, X; τ ) for the treatment confounding bridge function might be appropriate, where τ is an unknown finite dimensional parameter indexing the model, Corollary 2 suggests one can estimate τ by solving the estimating equation where m(W, A, X) is a user-specified function whose dimension is no smaller than τ 's.
Example 1 . If Z and W are both binary, rather than solving the system of equations implied by (12), one can instead specify a saturated model for the treatment confounding bridge function: and estimate τ = (τ 0 , τ 1 , τ 2 , τ 3 ) T by solving (14) with m(W, A) = (1, W, A, W A) T . Extension to Z and X with multiple categories is straightforward.
Example 2 . In case of continuous (U, X, Z), result (9) suggests the model If a univariate NCO W is available, we may solve (14) with m(W, A, X) = (1, W, A, X) T .

Estimation and Inference
In the previous sections, we have defined the causal parameter of interest β 0 as stratum-specific log risk ratio, introduced the treatment confounding bridge function as a key ingredient to identification of β 0 , and presented a strategy to estimate the treatment confounding bridge function leveraging an available NCO. We summarize the steps of our estimation framework in Algorithm 1 and present the large-sample properties of the resulting estimator ( β, τ ) in Theorem 3.
Algorithm 1 Negative control method to estimate vaccine effectiveness from a test-negative design 1: Identify the variables in the data according to Figure 1(e), in particular the NCEs and NCOs. 2: Estimate the treatment confounding bridge function by solving Equation (14) with a suitable parametric model q * (A, Z, X; τ ) and a user-specified function m(W, A, X). Write τ as the resulting estimate of τ . 3: Estimate β 0 by solving Equation (10) with q(A, Z, X) = q * (A, Z, X; τ ). The resulting estimator of β 0 is The estimated vaccine effectiveness is V E = 1 − exp( β).
Theorem 3 (Inference based on ( β, τ )). Under Assumptions 1-8 and suitable regularity conditions provided in Appendix J , the estimator ( β, τ ) in Algorithm 1, or equivalently, the solution to the estimating equation 1 n n i=1 G i (β, τ ) = 0 is regular and asymptotically linear with the i-th influence function The proof follows from standard estimating equation theory (See Van der Vaart (2000) Theorem 5.21). An immediate consequence of Theorem 3 is that we may estimate the variancecovariance matrix of ( β, τ ) with where HereÊ and V ar denote the expectation and variance with respect to the empirical distribution, respectively. A two-sided α-level Wald-type confidence interval of VE can then be obtained as where Σ n,1,1 is the (1, 1)-th entry of Σ n and z 1−α/2 is the (1 − α/2)-th quantile of a standard normal distribution.
The estimatorβ and the above confidence interval are constructed under the assumption that the disease is rare in the target population; for non-rare diseases,β is in general going to be biased and the confidence interval may not be well-calibrated. However, by Corollary 1, under the null hypothesis of no vaccine effects, the estimated q * (A, Z, X) converges to the true treatment confounding bridge function andβ is consistent for β 0 = 0. This implies that while our methods are approximately asymptotically unbiased for rare infections, they provide a valid test of no vaccine effect even if the infection is not rare.

Accounting for effect modification by measured confounders
So far we have operated under Assumption 3 that VE is constant across levels of (U, X). As we now show, this assumption can be relaxed to allow for potential effect modification with respect to X without compromising identification. This extension is particularly important because empirical evidence has indeed suggested that flu vaccine effectiveness may vary across sex and age groups (Chambers et al., 2018); and similar effect heterogeneity is of key interest in case of COVID-19 (Fernández Villalobos et al., 2021).
Instead of Assumption 3, we consider a less stringent assumption: Assumption 9 (No effect modification by unmeasured confounders).
where β 0 (x) are g(u, x) are unknown functions of x and u, x respectively.
Under condition 1, Assumption 9 further implies that β 0 ( x], i.e. the conditional causal RR as a function of x. Similar to Theorem 1, we have: Theorem 4. Under Assumptions 1, 2, 4, 5 and 9, we have that The proof of Theorem 4 is identical to that of Theorem 1 with β 0 A replaced with β 0 (X)A.
Identification and estimation of the treatment confounding bridge function are also essentially identical to that of Corollary 2. Therefore, it is straightforward to extend Algorithm 1 to allow effect modification by measured confounders. We describe the algorithm and the large sample properties of the resulting estimator in Appendix K.

Estimating VE under treatment-induced selection into TND sample
Thus far, unbiasedness of the estimating function V 0 has crucially relied on Assumption 2 that A does not have a direct effect on S. In some settings, the assumption may be violated if an infected person who is vaccinated is on average more likely to present to the ER than an unvaccinated infected person with similar symptoms, so that treatment or vaccination-induced selection into the TND sample is said to be present. In such settings, both estimatorsβ and β(X) produced by Algorithms 1 and 2 may be severely biased because Assumption 2 may no longer be valid. Crucially, we note that this form of selection bias can be present even in context of a randomized trial in which vaccination/treatment is assigned completely at random, if the outcome is ascertained using a TND, for example in the cluster-randomized test-negative design studies of community-level dengue intervention effectiveness Anders et al. (2018), Dufault and Jewell (2020), Jewell et al. (2019), and Wang et al. (2022). In this Section, we provide sufficient conditions for identification under treatment-induced selection. In this vein, consider the following assumptions: for a = 0, 1.
That is, the risk ratio association between infection status and selection into the TND sample is independent of vaccination status. Furthermore, Furthermore, identification relies on the following modified definition of a treatment confounding bridge function: Assumption 5 . There exists a treatment confounding bridge functionq such that for a = 0, 1, E[q(a, Z, X)|A = a, U, X] = 1/P (A = a|U, X, Y = 0, S = 1) almost surely.
Note that if the infection is rare in the target population in the sense of Assumption 7, then the treatment confounding bridge function defined in Assumption 5 in Section 2.3.1 satisfies (20) approximately .
We now introduce the identification of the OR with the following theorem: In addition to Assumptions 6, this last assumption requires that neither Y nor S is a causal effect of W . Figure 1(f) illustrates a DAG that satisfies our assumptions regarding (Z, W ). As can be verified in the graph, Assumption 6 is needed to ensure that collider stratification bias induced by the path A → [S = 1] ← W upon conditioning on S = 1 is no longer present.
Identification of the functionq is given below: Theorem 2 . Under assumptions 4, 5 and 6 , for a = 0, 1 we have that E[q(a, Z, X)|A = a, W, X, Y = 0, S = 1] = 1/P (A = a|W, X, Y = 0, S = 1) We prove Theorem 2 in Appendix M. As a result of Theorem 2 , the parameters in the treatment confounding bridge function can be estimated by solving moment equation (14).
In summary, the above discussion suggests that one can continue to use Algorithm 1 to estimate VE in presence of treatment induced selection bias, albeit on the OR scale and under a modified set of negative control conditions. Algorithm 2 can similarly be justified under treatment-induced selection with assumptions in this section, except that β 0 in Assumption 3 is replaced by the conditional log RR β 0 (X).
As a side note, Assumption 2 automatically holds under Assumption 2, and hence the above results in this section also apply to the setting in previous sections that is illustrated in Figure 1(e). We present this statement in the following corollary.

Simulation Study
To assess the empirical performance of our proposed method, we consider two settings with different types of confounding and negative control variables, and perform corresponding simulation studies.
In the first setting, we consider no measured confounder, a binary unmeasured confounder U , a binary NCE Z and a binary NCO W . To trigger selection among subjects with Y = 0, we let D be a binary indicator of the presence of other flu like illnesses. The treatment confounding bridge function is thus given by (8). We assume the distribution of Y is Bernoulli with a log- We consider values of β 0 to be -1.609, -0.693, -0.357 or 0, corresponding to a risk ratio of 0.2, 0.5, 0.7 or 1. We assume the selection S only equals one with nonzero probability if at least one of Y , W and D equals one, and is independent of A and Z conditional on other variables. The resulting prevalence of Y in the target population is 0.75% among the unvaccinated individuals and 0.55%, 0.65%, 0.72% or 0.75% among the vaccinated individuals, corresponding to four values of β 0 . Next, we consider a setting where X, U , Z and W are all univariate continuous variables. We generate the infection outcome using a log-linear model We generate A and Z following Example 2 in Section 2.3. As such the treatment confounding bridge function is given by Equation (9). The probability of S = 1 is 1 only if at least one of Y and D is nonzero. The resulting prevalence of Y in the target population is 0.34% among the unvaccinated individuals and 0.24%, 0.28%, 0.31% or 0.34% among the vaccinated individuals, corresponding to four values of β 0 . Appendix N and O give more details on the data-generating mechanism for the two settings.
In each scenario, we simulate a target population of size N = 7, 000, 000 and implement 1, 000 simulation iterations. For both settings, we evaluate the performance of three estimators for β 0 : • NC estimator: our proposed estimator given by Algorithm 1. In the first setting, we use a saturated parametric model (15) for the treatment confounding bridge function in the first setting, with m(W, A) = (1, W, A, W A) T ; in the second setting, we use model (16) and m(W, A, X) = (1, W, A, X) T .
• NC-Oracle estimator: the estimated treatment confounding bridge function in Algorithm 1 is only an approximation under Assumption 7, whose bias may affect the estimation for β 0 , as derived in Appendix I. We therefore include NC-Oracle estimator that uses the true treatment confounding bridge function q(A, Z, X). Appendices E include derivation of the true treatment confounding bridge function under the continuous (X, U, Z, W ) setting.
• Logistic regression: we also consider a logistic regression model of Y on A (and X in the second setting), overlooking the unmeasured confounders U . This is a common choice for covariate adjusted analyses of test-negative designs but ignores biases caused by U (Bond, Sullivan, and Cowling, 2016). We comment in Appendix P that the estimator is appropriate in the absence of unmeasured confounders except for potential model misspecification.
We note that the target parameter β 0 for NC estimator and NC-Oracle estimator is log causal RR, while logistic regression gives log causal OR. However, the two parameters are approximately equal under Assumption 7 where the infection risk is low in the target population. Figure 2 shows the bias of three estimators considered and the coverage of their 95% confidence intervals. In both settings, both NC and NC-Oracle are essentially unbiased whereas logistic regression gives a biased estimate in all scenarios. NC-Oracle exhibits slightly higher precision than NC, which implies that estimating the treatment confounding bridge function in the TND is only slightly more variable. The 95% confidence intervals for NC and NC-Oracle both achieve nominal coverage, whereas logistic regression-based confidence intervals under-cover severely. We repeated the simulation under a non-rare disease setting in Appendix Q. In such scenarios, while NC-Oracle estimator is still unbiased with calibrated 95% confidence intervals, the NC estimator is biased in general except when β 0 = 0. We conclude that the proposed NC estimator is unbiased of the log causal RR either under a rare disease setting or under a non-rare disease setting with no vaccine effect.

Application
We applied our proposed method to a TND study of COVID-19 VE against COVID-19 infec- We took immunization visits before December 2020 as NCE since COVID-19 vaccines were not available before December 2020 and immunization before was unlikely to affect the risk of COVID-19 infection; nor that of the selected NCOs we describe next. For NCO, we selected a binary indicator of having at least one of the following "negative control outcome" conditions after April 5, 2021: arm/leg cellulitis, eye/ear disorder, gastro-esophageal disease, atopic dermatitis, and injuries. Such candidate NCE and NCO are likely to satisfy the requisite conditional independence conditions for them to be valid negative control variables and to be related to a patient's latent HSB. We adjusted for age groups (<18, between 18 and 60, or ≥ 60), gender, race (white or non-white), Charlson comorbidity score ≥ 3, and the calendar month of a testpositive subject's first positive COVID test or a test-negative subject's last COVID test. Table 3 in Appendix R summarizes the distribution of negative control variables, demographic variables and COVID-19 infection among vaccinated and unvaccinated subjects.
Because NCE is expected not to be associated with either the outcome or NCE in a fully adjusted analysis unless there is unmeasured confounding, we first fit regression models to detect presence of residual confounding bias. Conditioning on the baseline covariates, in both vaccinated and unvaccinated groups, NCE is significantly associated with COVID-19 infection (p < 0.001) and NCO (p < 0.001) in corresponding adjusted logistic regression models, suggesting the presence of hidden biases (See Appendix R  (Hausman, 1978) comparing the two estimates on the log-RR scale is 8.84, giving a p-value of 0.003, indicating that the double negative control VE estimate is significantly smaller than that given by logistic regression.
The VE estimated with a standard logistic regression model is larger than VE estimates from RCTs which reported VE=94% for Moderna and VE=95% for Pfizer vaccines, respec-tively (Baden et al., 2021;Polack et al., 2020;Sadoff et al., 2021). We hypothesize that this may be due to some degree of confounding by HSB and related factors, which our proposed double NC approach appears to account for to some extent, recovering VE estimates more consistent with those of RCTs.

Discussion
In this article, we have introduced a statistical method for estimating vaccine effectiveness in a test-negative design. The approach leverages negative control variables to account for hidden bias due to residual confounding and/or selection mechanism into the TND sample. Negative control variables abound in practice, such as vaccination history which is routinely collected in insurance claims and electronic health records. Hence the proposed method may be particularly useful in such real world settings to obtain improved estimates of vaccine effectiveness.
The TND is a challenging setting in causal inference where selection bias and unmeasured confounding co-exist, selection is outcome-dependent, and unmeasured confounders also impact selection. As a result, the causal effect of interest is in general not identified from such studies (Cai and Kuroki, 2012). Nevertheless, we establish that progress can be made under a semiparametric multiplicative model, provided the outcome is rare in the target population, and double negative control variables are available. To this end, this article showcases the potential power of negative control methods and proximal causal inference in epidemiologic research (Shi, Miao, and Tchetgen Tchetgen, 2020;Tchetgen Tchetgen et al., 2020).
We focused on the outpatient TND, where recruitment is restricted to subjects who seek care voluntarily. TNDs have also been applied to inpatient settings for studying VE against, for example, flu hospitalization (Feng, Cowling, and Sullivan, 2016;Foppa et al., 2016). In inpatient TNDs, differential access to healthcare and underlying health characteristics between vaccinated and unvaccinated subjects are likely the main causes of confounding bias (Feng, Cowling, and Sullivan, 2016). Our methods are still applicable in such settings, but negative control variables should be selected to be relevant to the source of unmeasured confounding mechanism. For example, previous vaccination and hospitalization outside the flu season or hospitalization due to other flu-like illnesses are viable candidate NCE and NCO, respectively (Jackson et al., 2006).
Our approach is suitable for post-market TND studies where real-world vaccine effectiveness is of interest and vaccination history is obtained retrospectively, possibly through electronic health records. For vaccine efficacy in a controlled trial setting, Wang et al. (2022) recently developed estimation and inference of RR in cluster-randomized TND, aiming to correct for bias due to differential HSB induced by the intervention being unblinded. Because of randomization, they considered HSB as a post-treatment variable and proposed a log-contrast estimator which corrects for selection bias by leveraging a valid test-negative outcome, under an assumption that either (i) the vaccine does not have a causal effect in the population, and the causal impact of vaccination on selection is equal for test-positive and -negative subsamples; or (ii) among care seekers, the incidence of test-negative outcomes does not differ between vaccinated and unvaccinated, and the intervention effect among care seekers is generalizable to the whole population.
We note that even under randomization, identification conditions given in Section 2.6 are neither stronger nor weaker than those of Wang et al. (2022) described above, as neither set of assumptions appear to imply the other. An important advantage of our proposed methods is that they can be used to account for selection bias in a TND study irrespective of randomization.
Our methods target RR as a measure of VE instead of the more common OR (Jackson et al., 2006;Sullivan, Tchetgen Tchetgen, and Cowling, 2016). These two measures are approximately equal for rare infections. Schnitzer (2022) recently considered estimation of a marginal causal RR in the TND sample and justified the use of an inverse probability of treatment weighted (IPTW) estimator in a setting in which an unmeasured common cause of infection and selection into the TND sample does not cause vaccination (and thus there is no unmeasured confounding). Instead, our methods allow for an unmeasured common cause of vaccination, infection and selection into the TND sample; however in order to estimate a causal RR, we invoke both, an assumptions of no effect modification by an unmeasured confounder, and a rare-disease condition. As we establish, the latter assumption is not needed if there is no vaccine effect against infection outcome. In Section 2.6, we establish that under a homogeneous OR vaccine effect measure condition, and an alternative definition of the treatment bridge function, our methods can identify a causal effect of the vaccine on the odds ratio scale without invoking the rare disease condition.
Throughout the article, we have assumed diagnostic tests are accurate and individuals who seek care are sparsely distributed, such that the vaccination of a given subject in the TND sample does not protect another study subject from infection, , i.e. there is no interference in the TND sample, a common assumption in TND literature. This assumption may be violated if members of the same households present in the ER in which case block interference must be accounted using results from interference literature (Hudgens and Halloran, 2008;Tchetgen Tchetgen and VanderWeele, 2012). Sensitivity analysis may be considered to evaluate how violation of these assumptions can potentially bias inferences about VE.

A Proof of Proposition 1
We need to prove the result immediately follows if we can show that By Assumption 2, the left-hand side of (21) equals We further have =0.

B Existence of solutions to (7)
In this section, we provide the conditions of existence of solutions to (7). The conditions for (12) can be similarly derived. The results in this section directly adapted from Appendix B of Cui et al. (2020).

C Proof of Theorem 1
We need to prove it suffices to show that By Assumption 2, the left-hand side is The last equality is proved in Appendix A.

D Treatmeng bridge function and estimation with categorical NCE, NCO and unmeasured confounders.
We consider a categorical unmeasured confounder U with categories u 1 , . . . , u J and NCE Z with categories z 1 , . . . , z I . We assume there are no other covariates X. We write p ia.k = P (Z = z i , A = a|U = u j ) for i = 1, . . . , I, a = 0, 1, and j = 1, . . . , J.
The probabilities p ia.k 's can all be estimated from the study data.

E Derivation of the treatment confounding bridge function in Example 2
By Assumption 5, the treatment confounding bridge function q(A, Z, X) should satisfy for all a, u and X.
We write q(A, Z, X) = 1 + r(Z, A, X), then E[r(Z, a, X)|U, A = a, X] = 1 − P (A = a|U, X) P (A = a|U, X) . Consider This requires that we conclude the parameters in the bridge function are

F Proof of Theorem 2 and Corollary 1
We first introduce a few properties due to the rare disease assumption: Lemma 1. Under Assumptions 2, 4, 6 and 7, for every a, w, u and x, we have The rest follows similarly.
(b) For every a, u and x, we have By Lemma 1(a), we have The result follows by noticing that P (A = a|U = u, X = x, Y = 0) = P (A = a|U = u, X = x, Y = 0, S = 1) due to Assumption 2.

G Discussion of Assumption 8
Similar to (7), Equation (12) defines a Fredholm integral equation of the first kind, yet only involves observed data. Although a treatment confounding bridge function q(A, Z, X) must satisfy (12), there is no guarantee that solving (12) gives a treatment confounding bridge function if multiple solutions exist. When a solution to (12) exists, we give the following assumptions that guarantee the uniqueness of solution.
Assumption 10 (Completeness). were originally introduced by Lehmann and Scheffe to identify the so-called unbiase minimum risk estimator (Lehmann and Scheffé, 2012a,b). In econometrics and causal inference literature, completeness conditions have been employed to achieve identifiability for a variety of nonparametric or semiparametric models, such as instrumental variable regression (D'Haultfoeuille, 2011;Newey and Powell, 2003), measurement error models (Hu and Schennach, 2008), synthetic control (Shi et al., 2021), and previous works in negative control methods (Cui et al., 2020;Miao, Shi, and Tchetgen Tchetgen, 2018;Ying et al., 2021). The completeness condition holds for a wide range of distributions. Newey and Powell (2003) and D'Haultfoeuille (2011) provided justifications in exponential families and discrete distributions with finite support. Andrews (2011) constructed a broad class of bivariate distributions that satisfy the completeness condition.
(b) Under Assumptions 10(a) and 1-5, the treatment confounding bridge function is unique.
That is, if two square-integrable functions q(A, Z, X) and q 1 (A, Z, X) satisfy for all a, u, x almost surely, then q(A, Z, X) = q 1 (A, Z, X) almost surely.
We prove Proposition 2 below. Proposition 2 states that the completeness conditions in Assumption 10 and the definitions of the negative control variables lead to a third completeness condition that only involves the study data. Proposition 2 states the uniqueness of the treatment confounding bridge function, q(A, Z, X), and the solution to (12), q * (A, Z, X). The function q * (A, Z, X) can therefore be identified from the study data alone and is a good approximation of q(A, Z, X) by Theorem 2.

H Proof of Corollary 2 and further discussion
We first prove equation ( In fact, one can show that any regular and asymptotically normal estimator of τ that satisfies (12) has influence function of the form for an arbitrary function m(W, A, X). Therefore, any regular and asymptotically normal estimator of τ corresponds to the solution of the estimating equation (14) for some function m(W, A, X).
To prove this result, we see that for any parametric submodel that satisfies (12) and is indexed by s such that the true distribution corresponds to s = 0, we have m(W, A, X) S = 1 = 0 and thus Rearranging the terms, we have I Discussion on the rare disease assumption 7 We first describe a crucial identity that links the log risk ratio β 0 and the treatment confounding bridge function q(A, Z, X). Let q * be the function that satisfies (12): .
We introduce an additional regularity condition: Assumption 11 (Uniform continuity). For any fixed positive square-integrable function g(U ) and a small positive number 0 < η < 1, there exists some 0 < γ = γ(g, η) < 0 such that for almost all a and u. Let Under mild regularity, the estimator β is regular and asymptotically linear for β * 0 , and therefore It suffices to study and similarly, Therefore, we have .
We conclude that

J Regularity conditions and proof of Theorem 3
We denote τ * 0 as the true value of τ such that q(A, Z, X; τ * 0 ) = q * (A, Z, X). We will give the regularity conditions and proof that ( β, τ ) is a regular and asymptotically linear estimator of (β * 0 , τ * 0 ). Here β * 0 and q * are the biased versions of β 0 and q defined in Appendix I, respectively, although the biases are negligible when the infection is rare. Following Appendix I, β is also a regular and asymptotically linear estimator of β 0 if sup a,w,u,x P (Y = 1|A = a, W = w, U = u, X = x) < δ n , Assumption 11 holds and A set of regularity conditions are R.1 The function τ → q(A, Z, X; τ ) is Lipschitz in a neighborhood of τ * 0 ; that is, for every τ 1 and τ 2 in a neighborhood of τ 0 and a measurable functionq(A, Z, X) with E[q(A, Z, X)] < ∞, The condition R.1 and the fact that β → exp(−βA) is Lipschitz in a neighborhood of β * 0 imply the function (β, τ ) → M i (β, τ ) is Lipschitz in a neighborhood of (β * 0 , τ * 0 ) for every i. The remaining proof follows Van der Vaart (2000) Theorem 5.21.

K Estimating conditional causal RR in the presence of effect modification by measured confounders
Algorithm 2 below describes a straightforward extension of Algorithm 1 to estimation the conditional vaccine effectiveness V E(x) = 1 − exp(β 0 (x)) under Assumption 9 and a parametric model β 0 (X; α) indexed by a finite-dimensional parameter α.
Algorithm 2 Negative control method to estimate conditional vaccine effectiveness from a test-negative design 1: Identify the variables in the data according to Figure 1(c), in particular the NCEs and NCOs. 2: Estimate the treatment confounding bridge function by solving the equation (14) with a suitable parametric model q * (A, Z, X; τ ) and a user-specified function m(W, A, X). Write τ as the resulting estimate of τ . 3: Estimate α by solving Denote the resulting estimate of α as α. The estimated conditional vaccine effectiveness is We describe the large-sample properties of the estimator ( α, τ ) in the theorem below.
Theorem 5 (Inference of ( α, τ )). Under Assumptions 1-2, 4-5, 6, 9 and suitable regularity conditions listed at the end of this section, the estimator ( α, τ ) in Algorithm 2, or equivalently, the solution to the estimating equation 1 n n i=1G i (α, τ ) = 0 is regular and asymptotically linear with influence function Here C(X) is a user-specified function of X with the sample dimension as α.
Suppose that in Algorithm 2, one specifies β 0 (X; α) = X T α, then a natural choice for C(X) is C(X) = X. A sandwich estimator of the asymptotic variance of ( α, τ ) can be deduced from previous derivations. Under Assumption 9, we have shown that one can identify V E(X), however, one may be unable to identify the population marginal V E without an additional assumption. Interestingly, we note that the population marginal risk ratio would remain nonidentified even if one had access to a random sample from the target population to inform the marginal distribution of X. Specifically, as shown in Huitfeldt, Stensrud, and Suzuki (2019) Denote α * 0 as the value of α such that β 0 (x; α * 0 ) = β * 0 (x). Below we give the set of regularity conditions such that ( α, τ ) is a regular and asymptotically linear estimator of (α * 0 , τ * 0 ): R'.1 The function τ → q(A, Z, X; τ ) is Lipschitz in a neighborhood of τ * 0 and α → β 0 (X; α) is Lipschitz in a neighborhood of α * 0 . where c = a * ,y * P (A = a * |Y = 0, U, W, X, S = 1)P (Y = y * |A = 0, U, W, X, S = 1) exp(β 0 a * y * ).
Proof. Note that

N Simulation setting with binary unmeasured confounder
We generate the data of a size N = 7, 000, 000 general population according to Figure 1, with the following distribution: In the above data generating process, to mimic a test-negative design platform, we created a binary indicator D for flu-like diseases other than W . The study sample contains subjects with S = 1. The distribution of S indicates that only subjects with at least one of Y , D and W equal to one will be recruited into the study sample. We chose the parameters as in Table 1, which resulted in an average study sample size of between around 48,000 and around 52,000. Table 1: Parameter values of the data generating distribution in the simulation with a binary unmeasured confounder.

Parameter
Value To obtain the true treatment confounding bridge function, by (7), we have z q(a, z)f (z|u, a) = z q(a, z)f (a|u, z)f (z|u)/f (a|u) = 1/f (a|u) and thus for each u, a. We obtain that   q(0, 0)
We chose the parameter values according to Table 2, which results in an average study sample size of between around 43,000 and around 47,000. P Logistic regression as a naïve approximate estimator of log risk ratio, ignoring unmeasured confounders Figure 3 shows a causal diagram of a test-negative design with no unmeasured confounders.
Again we assumed selection into the study S is independent of the subjects' treatment status A, given the subjects' infection status Y and other covariates X.
A Y S X Figure 3: Directed acylic graph of a test-negative design with no unmeasured confounders.
In this scenario, we note that the conditional odds ratio given X in the study population equals the conditional odds ratio in the general population, i.e. This implies that we may estimate the stratum-specific odds ratio in the general population by fitting a logistic regression model to the study data: P (Y = 1|A, X, S = 1) = expit(γ 0 + γ X X + γ A A) where expit(x) = exp(x)/[1 + exp(x)], and γ A is the log odds ratio. When the outcome is rare in the general population, i.e. P (Y = 0|A = a, X = x) ≈ 1 for any a, x, the log odds ratio γ A also approximates the log risk ratio.

Q Simulation for non-rare diseases
To investigate the performance of our method for non-rare diseases, we repeat the simulation with the same setup with binary or continuous confounders. For the simulation with binary confounders, we set η 0Y to be log(0.20); for the simulation with continuous confounders, we set µ 0Y to be log(0.20) -both correspond to an infection risk of 20% for subjects with A = U = X = 0.
The results are in Figure 4. While the NC-Oracle estimator remains unbiased and maintains calibrated confidence intervals, the NC estimator is in general biased with under-covered confidence intervals. Notably, the NC estimator is unbiased with calibrated confidence intervals under the null hypothesis where β 0 = 0. The logistic regression estimator is biased and has unpredictable behavior.